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Graphene subjected to chiral-symmetric disorder is believed to host zero energy modes (ZEMs) 
resilient to localization, as snggested by the renormalization group analysis of the underlying non¬ 
linear sigma model. We report accurate quantum transport calculations in honeycomb lattices with 
in excess of 10® sites and fine meV resolutions. The Kubo dc conductivity of ZEMs induced by 
vacancy defects (chiral BDI class) is found to match 4e®/7rh within 1% accuracy, over a paramet¬ 
rically wide window of energy level broadenings and vacancy concentrations. Our results disclose 
an unprecedentedly robust metallic regime in graphene, providing strong evidence that the early 
held-theoretical picture for the BDI class is valid well beyond its controlled weak-coupling regime. 
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After more than half a century, Anderson localization 
remains a central concept in condensed matter physics, 
with its many ramifications providing new insights into 
the behavior of disordered electrons [1]. The discov¬ 
ery of the “tenfold” symmetry classes of disordered met¬ 
als [2, 3] —beyond the standard threefold Wigner-Dyson 
classification scheme—has revealed a surprisingly rich 
diagram of Anderson localization transitions, including 
multifractality and critical delocalization in low dimen¬ 
sions [4]. 

The interest in critical quantum transport in two- 
dimensional (2D) systems has been greatly amplified 
with the discovery of graphene, a one-atom-thick crys¬ 
tal endowed with massless Dirac fermions [5]. The in¬ 
ternal pseudospin of the Dirac fermions—stemming from 
the honeycomb lattice structure with two sublattices— 
enables a rich variety of quantum transport phenomena 
[6, 7], including minimum conductivity in the clean limit 
[8], and crossover from weak-localization—orthogonal 
class—to weak-antilocalization—symplectic class—with 
increasing impurity potential range [9]. 

Recently, disordered graphene in the chiral symme¬ 
try class has been the focus of much attention [10-13]. 
In chiral models defined on bipartite lattices, disordered 
wave functions come in electron-hole pairs with energies 
±E linked by a unitary matrix diagonal in the sublattice 
space, i.e., \<j)±) = az\4>^)- A remarkable feature of the 
chiral class is the existence of critical states at the band 
center—zero-energy modes (ZEMs)—possessing multi¬ 
fractal statistics and absence of weak localization correc¬ 
tions at all orders in perturbation theory [2]. In graphene, 
the simplest realization of critical ZEMs is provided by 
randomly distributed vacancies. A vacancy is a topo¬ 
logical defect obtained by cutting out all adjacent bonds 
to a given carbon site. Vacancies drastically affect the 
spectrum near the Dirac point, leading to the appear¬ 
ance of ZEMs with enhanced density of states (DOS) and 
quasilocalized character [14, 15], which can be detected 


by scanning tunneling microscopy [16]. Other examples 
of chiral-symmetric disorder in graphene include random 
non-Abelian gauge fields (ripples) [17], and resonant scat- 
terers (e.g., adsorbed hydrogen) [18]. Whether quantum 
criticality induced by chiral disorder could explain the 
resilience of the minimum conductivity of graphene to 
Anderson localization is an outstanding question. 

The focns of this Letter is on vacancy-induced ZEMs, 
recently implicated in a controversy regarding the exact 
nature of the quantum transport at the Dirac point [19- 
22] . Vacancy-defective graphene belongs to the chiral or¬ 
thogonal ensemble (class BDI in the Altland-Zirnbaner 
classification of random fermion models [3]). The van¬ 
ishing of the /3-function of the effective nonlinear sigma 
model (NLcrM) led Ostrovsky et al. to conjecture a line 
of fixed points with nonuniversal metallic conductivity of 
the order of the conductance quantum (t(0) « jh [22- 
24] . However, the validity of the NLuM of the BDI class 
has been questioned, as vacancies are infinitely strong 
scatterers, not amenable to perturbative analysis [12]. 
On the other hand, numerical evaluations of the con¬ 
ductivity using wave-packet propagation methods show 
localization of all states (j{E) 0, including the ZEMs 

[19-21]. The Gade singularity in the DOS approaching 
A —>• 0 [12], however, raises questions on the validity 
of the extraction of the conductivity using wave-packet 
propagation methods. 

In this Letter we report on accurate calculations of the 
longitudinal dc conductivity in macroscopic large disor¬ 
dered graphene. By employing an exact representation 
of the Kubo formula in terms of Chebyshev polynomials, 
we were able to extract the behavior of (7(E) at the Dirac 
point with unprecedented resolution. Our results univo- 
cally show that vacancy-induced ZEMs display critical 
delocalization, as suggested by perturbative calculations 
based on the NLcrM [22-24] and numerical studies of the 
two-terminal conductance in nanoribbons with resonant 
scalar impurities [22, 23]. We find a constant conductiv- 
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ity over a wide range of vacancy concentrations, 

4e^ 

cr(0) = CTZEM (1.00 ± 0.01) , CTZEM = ^ ■ 

Strikingly, the ZEM conductivity is found to be robust 
with respect to variations in the inelastic broadening 
parameter 77 entering in the disordered Green functions 
down to 77 = 2.5 meV. This result is very surprising as 
vacancies are the ultimate case of a strong short-range 
disorder in graphene mixing K and K' valleys [ 6 , 7]. 

The model .—Chiral disordered graphene is modeled by 
the standard tight-binding Hamiltonian of tt electrons de¬ 
fined on a honeycomb lattice 

H ={a]bj+b]a^'^ , ( 1 ) 

(*j> 



FIG. 1: Density of states of disordered graphene as function 
of Fermi energy. The Gade singularity of ZEMs is apparent as 
the energy levels are probed with increasing resolution 77 —>■ 0 . 
The pristine DOS is shown (black line) as a guide to the eye. 


where {i,j) denotes nearest-neighbor pairs of carbon 
atoms and t = 2.7 eV is the corresponding hopping in¬ 
tegral [5]. Periodic boundary conditions along zigzag 
and armchair directions are employed. The vacancies— 
obtained by removing the corresponding orbitals—are 
distributed randomly on both sublattices with overall 
concentration rii. In what follows, we briefly outline the 
Chebyshev-polynomial Green function method (CPGF) 
used to accurately evaluate spectral properties and re¬ 
sponse functions of real size systems. 

The CPGF approach. —The numerical evaluation of 
the lattice resolvent operator Q{z) = {z — H)~^ requires 
a nonzero broadening (resolution) parameter 77 = Im^: > 
SE, where 5E is the mean level spacing. We are inter¬ 
ested in the limit of small SE, where strong quantum 
interference effects associated with ZEMs can be fully 
appreciated [4]. Numerical evaluations of disordered lat¬ 
tice Green functions in the presence of critical states are 
computationally highly demanding. In Ref. [12] a time- 
domain stochastic method has been employed to extract 
the DOS with high resolution. Here, we evaluate target 
functions directly in the energy domain by expressing 
Green functions in terms of an exact polynomial expan¬ 
sion. Our approach turns out to be particularly advan¬ 
tageous in the calculation of the conductivity (see be¬ 
low). First-kind Ghebyshev polynomials {rn(a^)}raGNo are 
employed due to their superior convergence properties 
[25, 26]. The use of Ghebyshev polynomials as a basis 
set requires rescaling the spectrum of H into the interval 
[—1:1]. To this end, we scale both operators and energy 
variables, H ^ h = H/W, e = E/W, and A = 77 /IT, 
where W is the half-bandwidth. With this notation the 
Green function admits the following representation 

^ 00 

g{E + ip) = -Y.9n{c\)ru{^, (2) 

n—0 

where {Tn(h)} are dehned through the Ghebyshev recur¬ 
sion relations: 7o(h) = I, 7i(h) = h, and Tn+i{h) = 


2h -Tnih) — Tn-i{h). The coefficients {gn{e, A)}nGNo are 
system independent and possess a simple closed form [27]. 
The GPGF expansion (2) is the starting point of the ac¬ 
curate calculations reported in this work. 

Density of states .—We start with a brief discussion of 
the DOS. Formally, 

iy{E) = --^-TTlmg{E + ir]), (3) 


where Qs = 2 accounts for spin degeneracy and the bar 
means disorder averaging. According to Eqs. (2)-(3), the 
information about the DOS is contained in the Gheby¬ 
shev moments = Tr7/i(h) of individual disorder re¬ 
alizations. To probe features induced by chiral ZEMs 
with meV resolution, we consider a honeycomb lattice 
with D = 60 000 x 60 000 sites (wOd/rm^). This system 
has SE « 0.3 meV at the Dirac point in the absence of 
vacancies. The DOS for a dilute vacancy concentration 
Ui = 0.4% is shown in Fig. 1. Given the large size of 
the system simulated, one disorder configuration is suf¬ 
ficient to obtain very precise results. The expected en¬ 
hancement of the DOS associated with ZEMs near E = 0 
[14, 15] is seen to dramatically depend on the resolution. 
Extracting the exact scaling as E —)■ 0 is a demanding 
task as the number of Ghebyshev moments required to 
converge the DOS, i.e., N oc IT/ 77 , can be of the order of 
several tens of thousands even for meV resolution; here, 
TV = 15 X 10^. (Similar technical challenges were encoun¬ 
tered in Ref. [12].) The analysis of the data suggests that 
the singularity is stronger than that predicted by Gade 
and Wegner [2] in full consistency with the detailed nu¬ 
merical study of Ref. [12] and the analytical results in 
Ref. [13]; see Supplemental Material for full details [27]. 

Conductivity .—The finite-size Kubo formula reads 




Im ^(E -I- 777 ) v\\lm.g{E + ip) 7 )|| 


(4) 


where 7 }|| = [f\\,H]/ih is the velocity operator (taken 
along the zigzag direction) and ft is the area. Here, the 
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FIG. 2: Fully converged Kubo dc conductivity for a 0.4% 
vacancy concentration as a function of Fermi energy at se¬ 
lected values of rj. The calculation required = 6.4 x 10^ 
Chebyshev moments. The inset shows a zoom of the peak at 
the Dirac point. Statistical fluctuations of the data are within 
~ 1 %. 


broadening rj mimics the effect of uncorrelated inelastic 
scattering processes, thus defining a time scale = h/rj 
for phase coherence in the system [32, 33]. 

The calculation of (t{E) follows identical steps as out¬ 
lined for the DOS. The presence of two Green func¬ 
tions in Eq. (4) requires a double polynomial expansion, 
rendering the calculation computationally extremely de¬ 
manding. Analogously to the kernel polynomial method 
[18, 25], the full spectral information is now contained 
in the Chebyshev moments (Jnm = Tr [77i(^)f’||7^(^)fl||]- 
The number of moments required (= N"^) depends on 
the desired resolution. Typically, « 10 x {W/vi) con¬ 
verges the conductivity to two decimal places. From the 
knowledge of {anm} the dc conductivity cf(E) is quickly 
reconstructed. See Ref. [27] for details. 

Full spectral results .—We first provide a bird’s-eye 
view of (t(E) before specializing to the case of ZEMs. For 
modest resolutions, ^ 10 meV, the physically meaning¬ 
ful limit crn_>.oo(E) is achievable in relatively small sys¬ 
tems with D « 10^. The fully converged dc conductivity 
for a dilute vacancy concentration rii = 0.4% is shown in 
Fig. 2. The behavior of aa^oo{E) with decreasing ry (i.e., 
increasing r^) provides direct information on the quan¬ 
tum transport regime [e.g., lim,^_>o crn_>oo(£-) = 0(> 0) 
in the insulating (metallic) phase] [33] . The limit —>■ oo 
is implicit hereafter. In an energy window ~ ±0.2 eV 
around E = 0—excluding the Dirac point itself— <j{E) 
decreases as ry is lowered, showing that localization ef¬ 
fects become increasingly more important as the ther¬ 
modynamic limit T] —> SE —> 0 is approached. The effect 
is notably stronger in the vicinity of the Dirac point, 
where strong localization (cr < e^/h) takes place already 


for 7 y wlO meV. This indicates that the a priori unknown 
simulated inelastic lengths Li = Li{E,Ti) are sufficiently 
large that charge carriers can effectively experience local¬ 
ization. In contrast, at energies \E\ > 0.2 eV an increase 
of <j{E) with increasing n is observed. This suggests that 
at such energies the simulated Li is not yet sufficiently 
large to observe localization effects. This interpretation 
is further confirmed below. At the Dirac point, on the 
other hand, (j{E) seems insensitive to the inelastic broad¬ 
ening parameter, matching ctzem with 1% precision in 
the entire range (see inset to Fig. 2). The anomalous 
robustness of the dc conductivity as A —)• 0 is highly sug¬ 
gestive of a quantum critical point, in agreement with 
field-theoretical predictions [24]. 

High resolution results. —To probe the extension of de- 
localization effects at the Dirac point, we devise a scheme 
to enable the computation of a{E) with meV resolution. 
First, we recursively construct the vectors 

- OO 

\F±{E)) = —'^lm[gr,{t,\)]dl\ip), (5) 

n—0 

where \g)) = Y^^=iXi\i) is a real random vector, O" = 
Tn{h)v\\, and O" = v\\Tn(hi). The random variables {xi} 
are uncorrelated and taken from a uniform distribution 
with {{xi}) = 0. The series is truncated at n < N when 
convergence to the desired precision is achieved. Finally, 
the Kubo dc conductivity is obtained from 

2he^ 

a^{E) = — {^_{E)\^+{E)), (6) 

by averaging with respect to both disorder and random 
vector realizations, i.e., a{E) = {{a^{E))) [27]. We note 
that for ZEMs, Eq. (5) acquires a particular simple form, 
|(^±(0)) = W“^^„Im[ 52 n( 0 , A)]d^"|(p). The advan¬ 
tage of Eqs. (5) and ( 6 ) is that they do not require cal¬ 
culation of individual Chebyshev moments {anm} (cost 
oc N^). In practice, this allows us to reach fine resolution 
(higher N) and also much larger systems containing up 
to a few billion lattice sites [34]. 

The high-resolution conductivity data across the vari¬ 
ous transport regimes identified earlier is given in Fig. 3. 
For convenience, we define an effective system size L* = 
HttufIp as the length of a pristine graphene system hav¬ 
ing 6e = p at the Dirac point. The largest simulation 
has L* ~ 0.7 /im, corresponding to a broadening of only 
2.5 meV. The state vectors in Eq. (5) were calculated 
with N = 12 000 Chebyshev iterations. The ZEM con¬ 
ductivity shows no sign of localization, being numerically 
very close to uzem = 4e^/( tt/i) through a parametrically 
wide range of inelastic broadenings in the range [2.5, 60] 
meV. This is to be contrasted with the behavior of a{E) 
away from the band center. For instance, at energies 
E = {50,100} meV there is a strong suppression towards 
cr —)■ 0 as L* increases. The localization is stronger in the 
neighborhood of the critical point at zero energy, with 
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FIG. 3: Fully converged Kubo dc conductivity for a 0.4% 
vacancy concentration as a function of L*//imp at selected 
energies. Here limp — 2.24 nm is the average distance between 
vacancies. A large honeycomb lattice with 3.6 x 10® sites was 
simulated to obtain good precision at large L*. Statistical 
fluctuations of the data are within ~ 1%. 


states with E = 50 meV localizing first than those hav¬ 
ing E = 100 meV. This behavior can also be inferred 
from Fig. 2, which shows that the tendency as 77 —>■ 0 
(L* —>■ 00 ) is for states to localize first in the vicinity 
of the ZEMs. In the inset to Fig. 3 the behavior for an 
energy far away from the Dirac point is shown. A tran¬ 
sition from ballistic to localized regime is observed as L* 
increases. Eventually, as L* —)■ 00 , all states with E ^ 0 
become localized. The latter is consistent with the be¬ 
havior expected for random fermions in the BDI class 
[1, 4]. Crucially, however, our accurate numerical treat¬ 
ment shows that the chiral symmetry at if = 0 protects 
ZEMs from localization up to T* « 1 /rm. This exotic 2D 
metallic regime had been predicted by the renormaliza¬ 
tion group (RG) analysis of the NLctM for the BDI class 
[24], although a fully nonperturbative calculation of the 
microscopic conductivity able to capture strong quantum 
interference effects at the Dirac point was lacking until 
now. 

Universal ZEM conductivity .—We finally investigate 
the robustness of the ZEMs metallic conductivity against 
changes in vacancy concentration. According to the per¬ 
turbative RG analysis for white-noise disorder in the BDI 
class, (t(0 ) should depend weakly on the disorder strength 
[24]. The actual picture for vacancies—being infinitely 
strong scatterers—is difficult to predict based solely on 
field-theoretical methods [12, 35]. The little sensitivity 
of a{0) to the effective length L* intuitively suggests a 
small dependence with the defect concentration too. In¬ 
terestingly, numerical results for transport across narrow 
graphene strips show cr( 0 ) « ctzem with weak depen¬ 
dence on Ui [23], demonstrating that, although evanes¬ 
cent modes are strongly affected by scattering from va¬ 
cancy defects, the large number of modes available (large 
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FIG. 4: Impact of vacancy concentration on bulk transport. 
Top panel: localization of states with E = 0.1 eV as a func¬ 
tion of Lt at various vacancy concentrations. Bottom panel: 
variation of a{E) with Ui at selected energies. 


DOS) counteracts perfectly to restore graphene’s clean 
ballistic conductivity [ 8 ] . To investigate the possibility of 
a disorder-induced universal metallic regime in graphene, 
we perform accurate Kubo calculations over 2 orders of 
magnitude in n^. We take a fine broadening t] = 2.5 meV 
so as to guarantee that L* is sufficiently large to capture 
any marked localization trend near the Dirac point. Our 
results are summarized in Fig. 4. Away from the band 
center the conductivity is strongly decaying with Ui as 
expected. For instance, at E = 0.1 eV—a typical Fermi 
energy in experiments—the conductivity swiftly enters in 
the strong localized regime already for dilute concentra¬ 
tions Ui « 0.2%. The dependence of (j{E) with L* is well 
fitted by an exponential law a cx see top panel. 

(The dependence of with the defect concentration is 
shown in the inset to the bottom panel.) However, at the 
band centre ZEMs show no signs of localization even be¬ 
yond the very dilute limit up to concentrations n = 1 %. 
For completeness we provide the results for E = 0.4 eV 
where transport is ballistic in the simulated range of L* 
up to n « 0.8% (see also Fig. 3). 

We briefly comment on previous wave-packet propa¬ 
gation calculations reporting on cr(0) —>■ 0 [19-21]. The 
strong singularity of the DOS at if = 0 makes the nu¬ 
merical extraction of the conductivity from the Einstein 
relation for diffusive transport a-{E) cx I'iE) D{E) very 
challenging. Additionally, the level broadening inserted 
as the inverse of the time cutoff in the wave packet propa- 
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gation may not be equivalent to the broadening employed 
in the finite-size Kubo formula [Eq. (4)]. Although com¬ 
putationally much more demanding, our approach has 
the advantage of assessing directly the microscopic con¬ 
ductivity with no further assumptions. 

In summary, we have demonstrated critical delocal¬ 
ization of zero energy modes in graphene by means of 
accurate numerical evaluations of the Kubo conductiv¬ 
ity in real size disordered systems containing billions of 
carbon atoms. Rather remarkably, the absence of local¬ 
ization in the BDI class at the Dirac point is consistent 
with nonlinear sigma model predictions [24] and numer¬ 
ical studies of the Dirac equation [22, 23], suggesting an 
unprecedentedly robust metallic state in two dimensions. 
We hope that our work further encourages the use of 
accurate large-scale polynomial methods in the study of 
Anderson localization transitions. 
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We provide details on the Chebyshev-polynomial 
Green function (CPGF) method as well as a thorough 
description of the accurate large-scale numerical calcula¬ 
tions presented in the main text. 


I. CHEBYSHEV-POLYNOMIAL GREEN 
EUNCTION (CPGF) METHOD 


At the heart of the CPGF method is the exact ex¬ 
pansion of the Green function G{z) = (z — h)~^ for a 
disordered lattice in terms of first-kind Chebyshev poly¬ 
nomials [1]. Below, we provide a short description of 
their main properties and a brief derivation of the CPGF 
expansion. 

We assume that the spectrum of h falls in the inter¬ 
val I = [—1 : 1] [2]. Accordingly, in what follows, z 
is a rescaled complex energy variable, z := e + iX with 
A > 0. Chebyshev polynomials {'r„(a;)}„gNo satisfy the 
recursion relations 


To{x) = 1, Ti{x) = X, Tn+i{x) = 2xTn{x) - Tn-i{x), 

( 1 ) 

such that Tn{x) = cos (n arccos x). They obey the or¬ 
thogonality relations 

J dxUj{x)T„{x)Tmix) = ^ Sn,m , (2) 


where uj{x) = l/(7r-\/l — x^), thus forming a complete set 
in the domain I. For a function /(x) and x € X one can 
write the expansion 


fix) 


= ^{x) XI 

n=0 


1 + <5n,0 


Tnix), 


(3) 


where = Jj. dx f (x) T^ix). Upon truncation of the 
expansion, the Chebyshev polynomials distribute errors 
uniformly, providing a superior polynomial expansion 
with uniform resolution 6x oc 1/N, where N is the high¬ 
est polynomial order used [3]. 

Let {cjn} and {Im)} be the eigenvalues and eigenvec¬ 
tors of the Hamiltonian h. In order to hnd an exact 
expansion of the lattice Green function, 


G{e + zA) — X 

m 


\m) (m| 

e + * A — Cm 


(4) 


in terms of Chebyshev polynomials, we make use of the 
identity [4] 


2i 


= TXa— Mz)Tnix), |x| < 1, (5) 

rr(, 1 + On,o 


n—0 


where Jniz) is the Bessel function of order n, to recast 
(4) as 


1 r°° 

g(e + iX) = - 

« Jo 


n—O 


^71,0 


( 6 ) 

where {7^(h)} are operators dehned by the matrix ver¬ 
sion of the Chebyshev recursion relations (1), that is. 


To{h) = Id , Ti{h) = h, Tn+iih) = 2h-Tn{h)-Tn-iih ), 

(7) 

with D denoting the Hilbert space dimension. The 
Laplace transform of the Bessel function has a well- 
known solution [5] 

1 / ! - \ 71 

J - sj . (8) 


Using this expression, after the analytic continuation s —>■ 
—iz and some straightforward algebra one obtains [1] 


g{e -I- iA) = E 9nie + iX)Tnih ), 

n—0 


2i (z - i\/l - 

- 1 I A - - 

1 + On,0 V 1 — -2 


(9) 

( 10 ) 


In what follows we show how to use the CPGF expansion 
(9) to compute spectral properties of large systems. 


H. APPLICATION: DOS AND LONGITUDINAL 
DC CONDUCTIVITY 

The thermodynamic density of states (DOS) is for¬ 
mally given by p[e) = limA->.o liniD_>.oo A) where 

^(e, A) =-— Tr Im^(e-I-iA). (11) 

ttD 

Here bar denotes disorder averaging. Using Eq. (9) we 
easily obtain 


^ C.XJ 

l^i^,^) = -:^'^^H9nie + iX)]fin, ( 12 ) 

n—0 

with Chebyshev moments given by 

Pn = Tr[r„(h)]. (13) 

Similarly to the kernel polynomial method (KPM) [6], 
the calculation of the DOS amounts to the determina¬ 
tion of the Chebyshev moments. This scheme is very 
convenient as {/r„} can be efficiently calculated even in 
very large systems with modest computational resources. 
Once the moments are determined, the smeared DOS in 
the entire parameter space (e, A) can be quickly retrieved 
from Eq. (12). 
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In a practical calculation the expansion is truncated 
so as to obtain an order-7V approximation to the target 
function, 

1 w-i 

A) =^ Im [g„(e + iA)]^„. (14) 

n=0 


For not too small A, extremely accurate approximations 
can be obtained for modest N. However, the extrac¬ 
tion of the thermodynamic limit almost invariably re¬ 
quires large N. The recursive calculation of Chebyshev 
moments {nn} explore the matrix relations in Eq. (7) 
and is numerically very stable. Furthermore, the Cheby¬ 
shev expansion has well defined resolution [through the 
broadening parameter appearing in the Green function 
(4)]. These are two substantial advantages of polynomial 
methods as compared to, e.g., Lanczos recursion [ 6 ]. No¬ 
tice that kernel coefficients are absent in Eq. (14); thus, 
this expansion is not equivalent to that obtained through 
the KPM [1]. 

The convergence rate of the exact expansion (14) de¬ 
pends crucially on the smoothness of the target function. 
As shown by the authors in Ref. [7] the presence of sharp 
resonances in the DOS requires a particularly large num¬ 
ber of moments. As a rule of thumb, the number of 
Chebyshev moments N determine the resolution SeN ac¬ 
cording to SeN ~ 1/A^ [ 6 ]. For instance, to probe fea¬ 
tures with small width 77 an accurate calculation requires 
SeN Tj and hence many Chebyshev moments before the 
expansion (14) converges [7, 28]. 

Next we discuss the application of the CPGF to the 
calculation of the dc conductivity. The starting point is 
the finite-size Kubo formula at zero temperature. 


( 11 T 


Im Q{e + iX) Vx Im G{e + iX) 


(15) 


where Vx = [x,h]/iH is the velocity operator and D is 
the area. Here, the broadening parameter A defines a 
time scale Ti oc 1/A for phase coherence in the system 
[ 8 ]. Using Eq. (9) we easily find 


2He 


2 N-l 


^a^(cA) = ^ E Im[ffn(£ + *A)]Im[5m(c + fA)] Vnm, 


n,771=0 


where 


Vnm^ = Tr 


Xx (h) Vx 


(16) 

(17) 


The evaluation of the dc conductivity is computation¬ 
ally more demanding than the DOS due to the presence 
of a double sum in Eq. (16). The number of moments is 
now , which can severely limit the resolutions and/or 
system size attainable. However, as shown in what fol¬ 
lows, this limitation can be overcome if enough memory 
exists to store the random vectors used for a stochastic 
evaluation of the moments. 


III. EFFICIENT CALCULATION OF 
CHEBYSHEV MOMENTS 

The complexity of the trace evaluation in Eqs. (13) and 
(17) is 0{D^). However, for very large sparse matrices, 
such as those appearing in effective tight-binding models, 
the full trace can be replaced by a stochastic average. For 
instance, for the DOS one can replace Eq. (13) by 

1 ^ 

~ ^'^{r\Tn{h)\r), (18) 

r—1 

where |r) = complex random vectors with 

coefficients satisfying ((^i)) = 0 and = Sij (real 

vectors may be used for spin rotational and time reversal 
symmetric Hamiltonians) [9]. The number of operations 
required to compute (18) is now 0{D x R). 

It is often assumed that the error in (18) has the very 
favorable scaling 0(l/\/RD) [ 6 ]. However, for very large 
n the matrix Tn{h) is no longer sparse and a larger R (or 
a larger system size D) is needed to obtain a stochastic 
trace evaluation (STE) with good precision. In practice, 
for very large systems, with D « 10®, we found that a 
single random vector is enough to obtain errors below 1 % 
for n up to ten thousand. Details are given below. 

We now overview the recursive method that allows us 
to efficiently calculate Chebyshev moments. For a gen¬ 
eral introduction the reader is referred to the review by 
A. Weisse et al. [ 6 ]. For concreteness, we describe the 
calculation of conductivity moments, i.e., 

1 « 

Vnm = rj^'^{r\VxTnih)VxTTn{h)\r) . (19) 

r—1 

(The DOS moments are computed with a similar 
scheme.) Suppose we start with a random vector |r). 
Then, using the recursion relations [Eq. (7)], we obtain 

%n+i{h)\r) = 2hTmCh)\r) - r^_i(h)|r), (20) 

which inspires us to write 

where 

lr}m = Tm(h)lr) . (22) 

In the above, |r)o = |r). In fact, the best way to proceed 
is to define a second, auxiliary truncated basis {jr)}, "c = 
1,.., R, with 


lr}=Vxlr}. (23) 

Then, we can apply the Chebyshev recursion to write 

|f)„+i = 2/i |r)„ - . (24) 
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The idea now is to implement a recursive calculation for 
each pair of random vectors, {|r)} and {|r)}, to gener¬ 
ate two sequences {Ir)o,|?')Af-i} and {|f)o,|f)Ar_i}, 
since they can be used to directly compute the stochastic 
averages 

Vnm(r) = n{r\v^\r)m (25) 

needed for the calculation of the conductivity moments 
(19), i.e., 

1 ^ 

=-^V„^(r). (26) 

r—1 

If large amounts of RAM are available, one can recur¬ 
sively compute {!?’)«} and {|r)„} for all n, m = 0,..., N— 
1, store them, and then evaluate the coefficients Vnm(?') 
for each r [no need to store {[r),^} and {|r),r„} for more 
than a given r at any time] [10]. 

We now show how to evaluate efficiently the matrix 
elements V„m(?’) = using a site representation 

for the random vector. Let 

i’n\xk,yk) = {Xk,yk\r)n, (j)l^\xk,yk) = {xk,yk\r)m, 

(27) 

with k = where {xk,yk) are lattice site coordi¬ 

nates. Then, 

D 

Vnm{r)= [(l)^^\xk,yk)]*'ip^nHxk',yk')>^ 

x {xk,yk\vx\xk',yk') ■ (28) 


defects, tg = t/W where t is the carbon-carbon hopping 
integral and W is half-bandwidth). Clearly, there is no 
need to store the entire D x D Hamiltonian matrix; a 
connectivity table with information about neighbor coor¬ 
dinates {ni{xk,yk), ■■■,'nz{xk,yk)}k suffices. This shows 
that the current scheme is just limited by the mem¬ 
ory required to store the amplitudes 
needed to compute the overlap V„m(r) for any two vec- 
tors {|r)„,|r)^}. 

The calculation can be made substantially more efh- 
cient if we are just interested in evaluating the conductiv¬ 
ity in a small rectangular parametric grid {{cp} x {A^}}, 
1 < P,q < Pmax.gmax- The Chebyshev moments Vnm 
contain more information than any such grid since they 
allow one to retrieve the complete spectral conductivity 
according to Eq. (16). Recall that A is only limited by 
the number of Chebyshev iterations, min A oc N~^, and 
hence can be made arbitrary small by increasing N. The 
conductivity for each point in the grid can be calculated 
efficiently using the single-energy algorithm outlined in 
the main text. The idea is to write the conductivity 
o’Ar(ep, Aq) for each pair {cp, Ag} [see Eq. (16)] as 

CTAr(ep,Ag) = ^^^(v5frHep,Ag)|(^(([^(ep,Ag)), (32) 

r—1 

where 

N-l 

|<7’+^(ep, Ag)) = Y + i\)\vx\rn) (33) 

n—0 

and 


We can write 

{xkjyk\vx\xk',yk') = i{xk' - Xk){xk,yk\Mxk',yk') ■ (29) 


N-l 

\P-\^pAq)) = Y Im[5n(ep + *Ag)]|r„). (34) 

n=0 


When only nearest-neighbor hopping is allowed, there is 
a substantial reduction in the number of terms required 
to compute the matrix element: 

D 

Vnm{r) = iYl4>mHxk,yk)]*Y'^^'^^^^’^ + Tx,yk -I- Ty)x 
k — 1 T 

X {xk,yk\h\xk-\-Tx,yk + Ty) , (30) 

where the number of lattice vectors r depends on the 
topology of the problem. The number of computational 
steps is thus precisely Dxz, which is much lower than D^. 
In most cases of interest, the Hamiltonian matrix element 
is just a constant hopping amplitude —<s, in which case 
we have 

D 

Vnmir) = -its Y\4'mHxk,yk)]*Y,^^'^-^i^k+Tx,yk+Ty) ■ 
k — 1 T 

(31) 

Notice that tg is a dimensionless hopping amplitude since, 
by construction, || ± h|| < 1 (for graphene with vacancy 


Equations (33)-(34) can now be computed iteratively 
with only a few vectors stored in memory (instead of 
2x N vectors). The substantial reduction in memory al¬ 
location has allowed us to treat very large tight-binding 
systems, in excess of a billion atoms {D = 3.6 x 10®), 
with high resolution; see next section. 

IV. PARTICULAR CASE: GRAPHENE WITH 
RANDOM VACANCIES 

In this section we provide the full numerical details of 
the calculations presented in the main text. 

A. The DOS 

The DOS of a macroscopic large honeycomb lattice 
{Nx = Ny = 60000; periodic boundary conditions) 
with dilute randomly distributed vacancies (concentra¬ 
tion m = 0.4%) has been calculated using the CPGF 



4 


method and numerical implementations as described 
above. The A^-order approximation to the DOS is given 

by 


VN{E,r]) 


1 Ini[5n(e + zA)] 

R ^ ^ ttDW 

n—0 r—1 


{rn\rn) , 


(35) 


where W = 3t is graphene’s half-bandwidth and |r„) as 
defined in Eq. (22), E = eW, and ry = XW. The initial 
random vector used in the Chebyshev recursion reads 
as |ro) = where {xi} are generated from a 

uniform distribution on the interval [—-s/S,’x/S]- In such 


a large Hilbert space {D = x Ny = 3.6 x 10®), self 
averaging guarantees that a single random vector R = 1 
and one disorder realization suffice to obtain accurate 
results even for fine resolutions, that is, large N. 

The accuracy of the stochastic evaluation of z^Ar(e, A) 
is illustrated with a few examples in Table I. The su¬ 
perior precision [better than 0 . 1 % for zero-energy modes 
(ZEMs) investigated in the main text] is a consequence of 
the size of the system simulated. We note that at larger 
values of the resolution parameter A (rj) the data preci¬ 
sion improves because convergence is achieved at smaller 
values of N (see below). 



ZEM 

E = 0.05 eV 

E = 0.10 eV 

E = 0.20 eV 


1.0782 

2.1507x10'’^ 

1.6777x10"^ 

1.11326x10“^ 

52 

1.0786 

2.1495x10“^ 

1.6764x10"^ 

1.11259x10“’^ 

53 

1.0784 

2.1501x10“^ 

1.6705x10"^ 

1.11310x10“’^ 

max SiWsi — v\lv 

« 0.02% 

« 0.03% 

« 0.30% 

« 0.35% 


TABLE I: Estimation of the data precision. DOS [#states/(atom-eV)] for three independent system realizations (disorder and 
initial random vector |r)), labeled Si, S 2 , and S 3 , at several energies for a resolution rj = XW of 1 meV. The relative maximum 
deviation from the average u is shown in the last row. The estimated accuracy is confirmed below using a different approach. 


We now assess the convergence of the Worder approx¬ 
imation. N must be sufficiently large such that vpj(e,X) 
is well converged (say to 1 % accuracy or better) for the 
smallest desired resolution A. In Fig. 1, we show the vari¬ 
ation of the DOS of ZEMs A) with N [see Eq. (35)]. 
The calculations highlight the need for many thousands 
of Chebyshev iterations when the spectrum is probed 
with fine resolutions (i.e., a few meV). Similar conclu¬ 
sions hold for other energies (not shown). For compari¬ 
son we show the KPM approximation to the DOS using 
a Lorentz kernel [7, 25]. Despite being accurate in the 
limit N —>■ 00 , the KPM convergence rate is manifestly 
poorer in this case. 

Having established the convergence and accuracy of 
the CPGF method in the case of graphene with vacancies, 
we show the fully converged DOS for 1 meV resolution 
in Fig. 2. A single system realization and random vector 
was employed. As an independent error estimator we 
use the electron-hole asymmetry degree, i.e., |i^oo(e) A) — 
Voo{.—f-, A)|/i^oo(e 5 A). The magnitude of the error and its 
dependencies with the Fermi energy are consistent with 
the earlier statistical analysis. 

In order to illustrate the divergent behavior of the 
DOS at E = 0 we show in the right inset (Fig. 2) 
a plot of Ev{E,rj) at several values of the resolution. 
According to the standard nonlinear sigma model pic¬ 
ture [11], the thermodynamic DOS behaves as v{E, 0) —>■ 
exp[—I In lEll”^/®] as \E\ —>■ 0, whereas Hafner and 



N 

FIG. 1: Convergence of the A—order approximation to the 
DOS of ZEMs at selected values of resolution (broadening) 
parameter 77 = XW with VF = 8.1 eV. A single realization 
of a disordered system with = Ny = 60000 and 0.4% va¬ 
cancy concentration has been considered. The limiting value 
^'jv->oo(e, A) has been estimated—with precision better than 
1%—from the value of 7zjv(e, A) at A = 15000. The KPM 
result is shown (dotted line) for comparison. 


co-workers observed a stronger singularity v{E, 0) —>■ 
|£’|“^| In with 2 > cc > 1 [ 12 ] in consistency with a 

recent prediction [13]. Our results indicate Eu^E, 77 ) —>• 0 
for rj down to 1 meV, which is consistent with the numer¬ 
ical analysis of Ref. [12]. A more detailed analysis would 
be needed to reveal the exact dependence as obtained in 
the CPGF. 
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Our results show that the accurate determination of 
the spectral properties of disordered graphene is highly 
demanding, especially near the Gade singularity where 
fine resolutions are needed to capture the correct behav¬ 
ior. Similar challenges were reported in Ref. [12] where a 
time-domain stochastic method was used to extract the 
DOS. Finally, we note that the calculations are not sen¬ 
sitive to the system dimension as long as the mean level 
spacing is the smallest energy scale Se < X. This makes 
the CPGF a convenient tool to extract the thermody¬ 
namic limit. In what follows, we show how the GPGF 
behaves for the calculation of the Kubo formula. 



FIG. 2: DOS of graphene with vacancy defects {rii = 0.4%) 
as function of Fermi energy (green dashed line). The reso¬ 
lution of the calculation is r; =1 meV. A logarithmic scale 
has been chosen to highlight the singular behavior of i/(e. A) 
as e —>■ 0. The solid black line shows the DOS of pristine 
graphene as a guide to the eye. The insets show the esti¬ 
mated error as function of Fermi energy (left) and a close 
look at the DOS singularity at i? = 0 (right). The energy 
grid contains 1000 points. 


B. dc conductivity 

Below we provide the numerical details of the transport 
calculations presented in the main text. In Sec. B1 we 
focus on the full-spectrum algorithm used to produce the 
a versus E curve in Fig. 2 (main text). Details of the 
high-resolution calculations with D = 3.6 x 10® and N up 
to 12000 [Figs. 3 and 4 (main text)] are given in Sec. B2. 

1. Full spectral results 

As discussed in Sec. II, the knowledge of individual 
Ghebyshev moments Vnm = Tr [vxTn{h)vx%n{h)] en¬ 
ables the full spectral determination of the dc conduc¬ 
tivity. However, an efficient numerical implementation 


requires enough memory to store 2x N vectors of dimen¬ 
sion D [see Eqs. (25)-(25)], which in practice limits the 
attainable D and/or TV. To boost the size of the simu¬ 
lations we implemented the Ghebyshev recursive method 
(Sec. Ill) in machines with large RAM. Having the se¬ 
quences {lr)o,..., |r’)Ar_i} and {lf)o,..., |f')Af-i} stored in 
RAM allows for a quick evaluation of the Ghebyshev mo¬ 
ments through optimized linear algebra subroutines. 



F(eV) 



N 


FIG. 3: Analysis of full-spectral results. Top panel. Gon- 
ductivity of graphene with vacancy defects {rii = 0.4%) as 
function of Fermi energy. The resolution of the calculation is 
p =10 meV [a is given in units (Jzem = 4e^/(7r/i)]. The en¬ 
ergy grids contain 1000 points. The inset shows the estimated 
error based on the standard deviation of 20 independent sets 
(each containing an average over 250 random vectors). For 
comparison, the error estimated using the electron-hole asym¬ 
metry degree is also shown. For clarity the grid in the inset 
contains only 49 points with E > 0. Bottom Panel. The 
convergence of the A-order approximation to the ZEMs mi¬ 
croscopic conductivity crjv(0,?7) is shown at selected values 
of 77. Glearly, several thousands Ghebyshev iterations (cor¬ 
responding to tens of millions expansion moments Vnm) are 
required as rj enters the meV range. 
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FIG. 4: Numerical analysis of the large-scale calculations. 
Top panel. “Single-shot” ZEM conductivity of realistic size 
graphene samples with a dilute concentration of vacancies 
(ui = 0.05%) as function of N. The resolution of the calcula¬ 
tion is r) =5 meV. Each data set (blue circles) corresponds to 
a single system realization. The average over 20 independent 
system realizations is shown in squares. Bottom Panel. The 
same as in the left panel for N = 12000 polynomials and a 
resolution rj =2.5 meV. 

The results reported in this section are for a hon¬ 
eycomb lattice with D = 3200 x 3200 sites and a va¬ 
cancy concentration = 0.4%. In order to extract the 
Kubo conductivity with satisfactory resolution we com¬ 
puted N = 8000 Chebyshev iterations [corresponding to 
iV^ = 6.4 X 10^ moments in the expansion of aN{E,ri), 
Eq. (16)]. The resulting N x N matrix is subsequently 
used to evaluate aN{E,r]) on a fine grid. D is large 
enough so that the thermodynamic limit 12 —> oo can 
be safely extrapolated. In Fig. 3 (a) we show <j{E) = 
limn_>oo <^N^oo(E,rj = lOmeV) for a fixed disorder real¬ 
ization. Here, E is the Fermi energy in eV. 

Remarkably, the stochastic trace in Eq. (19) required 
thousands random vectors to converge um^E, rf) to a good 
precision [14]. The high degree of electron-hole symme¬ 
try a{E,r]) = a{—E,ri) achieved [see Fig. 3 (a)] testifies 
to the high quality of the results. The error in aN{E,X) 
is estimated to be in the range 0.1-1%. This is further 
confirmed with a detailed numerical study summarized 
in the inset to Fig. 3 (see caption for details). 

The convergence of the 7V-order approximation for 


ZEMs is shown in Fig. 3 (b). Whereas for poor reso¬ 
lutions « 20 meV a few thousand Chebyshev iterations 
are sufficient, probing resolutions « 1 meV is manifestly 
more demanding. Moreover, statistical fluctuations in 
the STE become important at small 77 , which requires 
more random vector realizations (see the noise in the 
curve for 77 = 5 meV). Importantly, all the curves studied 
converge to ctzem to 1% accuracy, the main result of the 
Letter. A dedicated calculation at A = 0 will confirm 
this (see below). 


2. Single-energy high-resolution results 

In Sec. Ill we devised a “single-energy algorithm” that 
bypasses the computation of Chebyshev moments Vnm, 
allowing us to reach much larger system sizes. We now 
describe its application to the problem of the ZEMs in 
graphene. The calculations summarized in this section 
are for a honeycomb lattice with D = 3.6 x 10® sites. 
The huge system dimension results in a{E, 77 ) data with 
satisfactory accuracy even for a single system realization, 
i.e., one random vector i? = 1 and a single (vacancy) dis¬ 
order realization. This situation is computationally very 
convenient as it provides a quick “single-shot” evaluation 
of the dc conductivity. 

In Fig. 4 we show the variation of cr a 7 ( 0 , 77 ) with N. As 
mentioned, a single system realization converges fTAr( 0 , 77 ) 
to a very reasonable precision (note that the vertical axis 
is zoomed around a = ctzem)- The error bars increase 
slowly with N as the matrices Tn(h) become less and less 
sparse as n —)■ A^ — 1 for A^ 3> 1. For a dilute vacancy 
concentration rii = 0.05% and broadening 77 = 5 meV, 
we obtain (crAr= 8 OOo( 0 )) = 1.008 (in units of ctzem) with 
standard deviation 5a = 0.009 (corresponding to 0.8% 
of the mean value). For the high resolution calculations 
(77 = 2.5 meV), these values are (crAr=i 2 ooo( 0 )) = 1-006 
and 0.016, respectively. 

In the main text, a set of “single-shot” calculations 
with 77 = {2.5,5,7.5,10,12.5,15,20,40,60} meV and 
Ui = 0.4% (Fig. 3), and 77 = 5 meV and rii = 
{0.05,0.1,0.2,0.4,0.6,0.8,1}% (Fig. 4) were presented. 
In order obtain a conservative estimate of the error bars 
involved we performed 20 independent realizations of the 
more disordered system, i.e., rii = 1%- We obtained 
(crAr=8OOo(0)) = 1.014 with standard deviation 0.007, 
which suggests an accuracy of « 1 %. 

Probing resolutions resolutions approaching 1 meV be¬ 
comes increasingly more challenging as the number of 
iterations N increases considerably, and hence the num¬ 
ber of random vectors necessary to converge the STE. 
We performed a small set of simulations for 77 = 1 meV 
{N = 12000, averaged over 2 disorder realizations and 10 
random vectors) and obtained an average 0.95(Tzem with 
5% standard deviation. 
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